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1.  Introduction 


Many  of  the  inventory  models  which  are  used  in  practice  rely  upon  knowing 
the  probability  distribution  of  demand  over  a  leadtime.  The  common  assumption 
is  that  this  distribution  is  normal.  However,  in  certain  circumstances,  the 
normality  assumption  may  be  inappropriate.  The  purpose  of  this  paper  is  to 
derive  the  exact  distribution  of  leadtime  demand  under  the  following  assumptions: 
customer  requisitions  occur  according  to  a  stationary  Poisson  process,  requisition 
sizes  follow  a  logarithmic  distribution  and  leadtime  is  a  random  variable  with 
the  gamma  distribution.  In  addition  to  deriving  the  exact  distribution  of  lead- 
time  demand,  we  compare  our  results  to  actual  operational  data  and  discuss  a 
variety  of  approximations. 

A  number  of  researchers  have  considered  the  problem  of  determining  inventory 
operating  policies  when  requisition  size  exceeds  one.  For  example,  Hausman  [6] 
extends  Hadley  and  Whitin's  [5]  heuristic  while  Archibald  and  Silver  [1]  derive 
optimal  (s,  S)  policies.  These  studies  differ  from  ours  in  two  ways.  First,  in 
every  case  leadtirae  was  assumed  to  be  deterministic.  Second,  they  focus  primarily 
on  describing  optimal  and  suboptiraal  ordering  policies.  Our  interest  is  in  a  de¬ 
tailed  examination  of  the  distribution  of  demand  over  leadtime. 

2.  The  Logarithmic  Distribution 

The  logarithmic  (or  log  series)  distribution  was  originally  derived  by 
Fisher  et.  al.  [4]  and  has  been  discussed  by  Sherbrooke  [11]  in  connection  with 
inventory  problems.  It  can  be  derived  as  a  limiting  case  of  the  negative  bi¬ 
nomial  distribution  and  has  the  form 

(1)  f  (x)  *  0^ _  for  x  =  1,  2 

-x  ln(l-G) 


where  0  <  0  <  1.  Chakrauarti  et.  al.  [3]  recommend  the  method  of  moments  be 
used  to  estimate  0.  It  is  easy  to  show  that 

(2)  E(x)  =  0 _ 

-(1-0) In (1-0) 

A 

which  means  that  an  estimator  for  0,  say  0,  solves  the  trancendental  equation 

(3)  X  -  © _ 

-(1-0) In (1-0) 

where  x  is  the  observed  sample  mean.  Since  the  right  hand  side  is  an  increasing 
function  of  0,  this  equation  can  be  solved  very  efficiently  by  interval  bisection. 

We  have  collected  data  describing  the  requisition  size  distribution  for  a 
number  of  Air  Force  FOQ-type  (i.e.  consumable)  items.  For  many  of  these  items, 
the  logarithmic  distribution  appears  to  be  a  useful  approximation  to  the  observed 
data.  An  example  of  a  specific  item  is  presented  in  Table  1.  In  this  case, 

A 

the  observed  sample  mean  is  3.94,  which  results  in  0  =  .901.  Notice  the  very 
close  agreement  between  the  observed  and  the  predicted  cumulative  distribution 
functions  for  this  item. 

3.  The  LPG  Distribution 

Let  us  now  assume  that  requisitions  are  generated  by  a  Poisson  process  and 
the  requisition  size  has  a  logarithmic  distribution.  (That  is,  the  demand  pro¬ 
cess  is  a  compound  Poisson  process  with  logarithmic  compounding  distribution). 

It  is  well  known  that  the  total  number  of  units  demanded  in  any  fixed  time,  t, 
say  Z(t),  has  the  negative  binomial  distribution.  In  particular,  we  obtain 

(4)  P{Z(t)  *  x}  *  (ct  +  x  —  1)!  (1-Q)ct  0  for  x  *  0,  1,  2,  ... 

x!  (ct-1) ! 

where  c  *  -A/ln(l-0)  and  A  is  the  requisition  arrival  rate. 
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Table 

1.  Comparison  of 

observed  frequencies  and  those  predicted  by 

logarithmic 

distribution 

for  a  typical 

E0Q  type  item. 

X 

Number  of 
Observations 

Observed 

Frequency 

Theoretical 

Frequency 

Observed 

Cumulative 

Theoretical 

Cumulative 

1 

93 

.4247 

.3896 

.4247 

.3896 

2 

31 

.1416 

.1755 

.5663 

.5651 

3 

13 

10594 

.1054 

.6257 

.6705 

4 

15 

.0685 

.0712 

.6942 

.7417 

5 

10 

.0457 

.0514 

.7399 

.7931 

6 

15 

.0685 

.0386 

.8084 

.8317 

7 

8 

.0365 

.0298 

.8449 

.8615 

8 

8 

.0365 

.0235 

.8814 

.8850 

9 

3 

.0137 

.0188 

.8951 

.9031 

10 

4 

.0183 

.0152 

.9134 

.9190 

11 

7 

.0320 

.0125 

.9454 

.9315 

12 

3 

.0137 

.0103 

.9591 

.9418 

13 

0 

.0000 

.0086 

.9591 

.9504 

14 

1 

.0046 

.0072 

.9637 

.9576 

15 

2 

.  .0091 

.0061 

.9728 

.9637 

16 

1 

.0046 

.0052 

.9774 

.9689 

17 

0 

.0000 

.0043 

.9774 

.9732 

18 

1 

.0046 

.0037 

.9820 

.9769 

19 

0 

.0000 

.0031 

.9820 

.9800 

20 

• 

• 

• 

2 

• 

• 

• 

.0091 

• 

• 

.0027 

# 

• 

• 

.9911 

• 

• 

• 

.9827 

• 

• 

• 

2  .0091  .0008  1.0000  .9915 


This  result  appears  to  be  due  to  Quenouille  [10].  Baswell  and  Patil  [2] 
give  fifteen  different  derivations  of  the  negative  binomial  distribution,  thus 


accounting  for  its  power  in  describing  many  common  phenomena. 

Now  let  us  assume  that  the  procurement  leadtime,  x,  is  a  continuous  non¬ 
negative  random  variable  with  probability  density  g(x).  In  general,  the  number 
of  units  demanded  in  time  T  is  a  random  variable  with  probability  function  h(x) 
given  by 

oo  * 

(5)  h(x)  =  j  f  (x J  T )  g (t )  di 

o 


where  f(x|x)  is  the  probability  function  of  the  number  of  units  demanded  in  a 
time  x.  Under  our  assumptions,  f(xlx)  has  the  negative  binomial  distribution. 
Since  cX  is  in  general  not  an  integer ,  we  use  the  gamma  function  representation 
for  the  factorials,  so  that 


(6)  h(x) 


-£  7 

V*  ‘r. 


.  CT 


r(cx  +  x) (1-0)  g(T)  dT. 
'o  r(ct) 


Using  the  fact  that  T (a)  =  (a-1)  T  (a-1)  we  have 


x-1 

(7)  r(cx  +  x)  =  n  (ex  +  j) 

r(ct)  j-o 


x  k 
E(ct)kS 

k=l 


xk* 


where  the  coefficients  are  known  as  Stirling  numbers  of  the  first  kind  and 
can  be  computed  from  the  recursion 

<8>  s*k  ■  s»-l,k-l  +  <*-l>  sx-l,k’ 

for  k  3  ,  2,  •  •  • ,  x  and  x  *  1,  2,  •••  } 

with  S  -  0  for  all  x. 
xo 

Furthermore,  from  the  definition  of  c, 

(9)  (1-0) cT  ■  exp{cT  In  (1-0)}  -  e  ^T, 

so  that  we  may  now  write 
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(10) 


h(x)  =  0^ 

* 


X 

E 

k=l 


Ou 

f  k  -Ax  .  , 

J  x  e  g(x)  d i. 


We  now  specialize  to  the  case  where  g(x)  has  the  gamma  distribution  with 
parameters  ot  and  3  so  that 

(11)  g(x)  =  3  x  e  for  t  ^  o. 

r(a) 

Since  leadtimes  must  be  non-negative,  the  gamma  distribution  should  provide 
sufficient  flexibility  to  model  leadtime  variability  in  many  operating  environ¬ 
ments. 

Using  the  fact  that 

00 

(12)  /  Tk+ot_1  e-(A+3)f  =  p(k+q) 

o  k+a 

(A+6) 

and  that,  as  above, 
k 

(13)  I'(k+q)  =  E  aJS  , 

r<co  j=i  kJ 

we  obtain  the  following  as  the  probability  function  for  the  number  of  units 
demanded  in  a  leadtime: 

x  k. 

(14)  hu)  \a  0^  E  ( c  \  k  S  E  aj  S.  for  x  =  1,  2,  3,  ... 

\£+@ J  x!  k=l  \A+g/  j=l  J 

and  h(o)  ■/  a. 

V+0/ 

We  call  this  the  LPG  distribution  (for  Logarithmic-Poisson-Gamma) .  Its 
four  parameters  are  a,  6,  0  and  A  (0  and  A  determine  c).  An  example  of  the  LPG 
distribution  is  presented  in  Table  2  for  a  *  1,  3  *  1,  G  *  .8  and  A  *  1. 


4.  A  Recursion  for  Integer  a 
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For  numeric  calculations,  we  found  the  following  recursion  to  be  useful. 
Let  us  assume  that  a  is  integer,  and  let  =  (g/(A+g))u  and  C  =  c/(A+g). 
Further,  using  (13)  let  us  define 


(i5>  Tx,k  -  ci  •  K  c2  s=.,k 

x!  (a  -1)! 


Hence,  h(x)  may  be  computed  as  a  sum  of  T 


x,k- 


(16)  h(x>  -  Ek.1  Tx_t 


Note  that  T  =0  since  S  =0 
x,o  x,o 


and  that 


Tx,x  •  C1  •  ^  <<*  + 

x!  (a  -  1)! 


<17>  ■  fS  <“  +  *-»  Vl.x-1 

X 


since  S  =1.  Using  (8),  we  may  now  write  T  in  terms  of  T  _ 

x ,  x  x ,  K  x—r ,  k 

specifically,  we  obtain: 


(18) 


x,k 


0  (C2  (a  +  k  -  1)  T  J  ♦  (x  -  1)  Vl>k) 


Thus,  h(x)  may  be  evaluated  using  only  ^  terms.  This  provides 

significant  reductions  in  computer  memory  and  calculation  requirements  compared 
to  a  direct  evaluation  of  (14)  for  each  x. 


5.  Approximations 

Many  inventory  models  require  computing  reorder  points  from  fractiles  of 
the  leadtime  demand  distribution.  Finding  exact  fractiles  of  the  LPG  distri¬ 
bution  might  be  too  demanding  computationally  for  many  real  applications.  In 
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Table  2. 

The  L 'G 

Distribution 

Pararae  ters 

0-1,  B« 

1,  0-.8,  A=l. 

X 

h(x) 

H  (x) 

0 

.5000 

.5000 

1 

.1243 

.6243 

2 

.0806 

.7049 

3 

.0589 

.7638 

4 

.0451 

.8089 

3 

.0355 

.8444 

6 

.0283 

.8727 

7 

.0228 

.8955 

8 

.0J  85 

.9140 

9 

.0151 

.9291 

10 

.0124 

.9415 

11 

.0101 

.9516 

12 

.0083 

.9599 

13 

.0069 

.9668 

14 

.0057 

.9725 

15 

.0047 

.9772 

16 

.0039 

.9811 

17 

.0032 

.9843 

18 

.0027 

.9869 

19 

.0022 

.9892 

20 

.0018 

.9910 

i; 

i: 


to 
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this  section  we  consider  an  approximation  which  uses  a  scaled  version  of  the 
Poisson  distribution  to  approximate  the  negative  binomial  distribution. 

The  mean  and  variance  of  Z(t),  the  number  of  units  demanded  in  time  t, 
are  respectively 

(19)  E(Z(t>)  =  ct0/(l-0). 

(20)  VAR (Z  ( t )  )  =  ct0/(l~0)2, 

which  gives  VAR(Z\t) ) /E(Z(t ) )  =  ( 1—0)  ^  .  In  certain  circumstances,  one  may  have 
knowledge  of  the  variance  to  mean  ratio  of  the  demand  which  can  then  be  used  to 
estimate  0  directly. 

The  approximation  is  based  upon  replacing  the  negative  binomial  distribution 
of  Z(t)  with  a  scaled  Poisson  distribution.  Let  Y  be  a  Poisson  random  variable 
with  parameter  fit  and  let  W  be  defined  by  W  =  kY  for  some  k  >  o.  We  may  think 
of  W  as  a  random  variable  which  assumes  values  0*  k,  2k,  ...  and  whose  distri¬ 
bution  depends  upon  the  two  parameters  pt  and  k. 

Since 

(21)  E  (W)  =  kpt 

(22)  VAR(W)  =  k2pt 

we  have  VAR(W)/E(W)  =  k^  Thus,  we  set  k  =  (1  -  0)  ^  to  achieve  the  same  variance 
to  mean  ratio.  Comparing  the  mean  and  variances  of  W  and  Z(t)  we  see  that  p  =  cG 
(recall  that  c  =  -A/ln(l-0)). 

Since  the  negative  binomial  distribution  is  defined  on  all  non-negative 
integers,  we  would  like  the  approximation  to  be  defined  on  the  non-negative 
integers  as  well.  We  have  found  the  following  procedure  works  well.  Assume 
that  the  scaled  Poisson  probabilities  are  shifted  to  k/2,  3k/2,  5k/2,  ...  so 


that 
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(23)  PfW  =  (n+1 )  k/2 l  =  e~MMn  u=0,  1,  2,... 

n! 

We  then  assume  that  the  cumulative  distribution  function  is  linear  between  nk/2 
and  (n+l)k/2.  As  an  example,  suppose  0  =  .75,  t  =  1,  c  =  2  (that  is  A  =  2.77). 
Then  u  =  1 .5,  k  =  4  and 

P{w  =  21  =  e  ~Ut  =  .2231 

P{W  =  6 1  =  e_tJlnt  =  .  3347 

PlW  =  10}  =  e"llt(ut)2/2  =  .2510 

P(W  =  14}  =  =  .1255 

etc, 

The  comparison  of  the  exact  negative  binomial  probabilities  and  the  scaled 
Poisson  approximation  is  presented  in  Table  3  for  this  case. 

We  now  obtain  an  approximation  to  the  LPG  distribution  by  averaging  the 
scaled  Poisson  approximation  of  the  negative  binomial  with  the  gamma  distribution 
of  leadtime.  That  is, 

oo 

(24)  P  { Z  ( T )  =  x}  *  /  e~P'r(uT)x/k  BaT°t~1e~BT  di. 

o  (x/k)!  '  T(a) 

But  this  integral  is  exactly  a  Poisson  mixture  with  a  gamma  distribution 
which  is  still  another  way  that  the  negative  binomial  distribution  can  be  derived 
(see  Baswell  and  Patil  [2]).  Hence,  the  approximation  for  the  LPG  distribution 
is  a  scaled  version  of  the  negative  binomial  distribution.  The  approximation 


for  x  *  0,  1,  2,  . . . . 


therefore  is: 
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i 

i 


Table  3.  Comparison  of  Negative  Binomial  and  Scaled  Poisson  Approximations 

(0.75,  t=l,  c=2 ,  A=2 .77) 


X 

Negative  Binomial 
Probabil i ties 

Negative  Binomial 
Cumul.  Probabilities 

Scaled  Poisson 
Cumul.  Probabilities 

Scaled  Poisson 
Probabilities 

0 

.0625 

.0623 

.0744 

1 

.0938 

.1363 

.0744 

2 

.1055 

.2618 

.2231 

.0744 

3 

.1055 

.  3673 

.0837 

4 

.0989 

.4662 

.0837 

3 

.0890 

.3332 

.0837 

6 

.0779 

.6331 

.5578 

.0837 

7 

.0667 

.6998 

.0628 

8 

.0563 

.  7561 

.0628 

9 

.0469 

.8030 

.0628 

10 

.0387 

.8417 

.8088 

.0628 

11 

.0312 

.8729 

.0314 

12 

.0257 

.8986 

.0314 

13 

.0208 

.9134 

.0314 

14 

.0167 

.9361 

,9343 

.0314 

13 

.0134 

.9495 

.0118 

16 

.0106 

.9601 

.0118 

17 

.0085 

.9685 

.0118 

18 

.0667 

.9753 

.9814 

.0118 

19 

.0053 

•  9806 

.0035 

20 

.0040 

•  9846 

.0035 

21 

.0033 

.9879 

.0035 

22 

.0026 

.9905 

.9955 

.0035 

j 
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Note  that  these  probabi li t ies  are  defined  on  0,  k,  2k,  .  ...  As  with  the 
scaled  Poisson  we  suggest  shifting  these  probabilities  to  kl2,  3kl2,  ...  and 
approximating  the  probability  function  by  assuming  the  cumulative  distribution 
function  is  linear  between  these  f Tactile  points.  We  tested  a  variety  of  cases 
and  found  the  fit  to  be  excellent,  especially  in  the  tails.  In  Table  4  we  com¬ 
pare  the  exact  LPG  probabilities  for  the  parameter  set  considered  in  Table  2 
with  the  scaled  negative  binomial  approximation.  Note  that  since  0=.8,  we  have 
k-5  and  the  approximate  cumulative  probabil i ties  (labelled  H(x)  in  the  table) 
are  defined  at  the  points  2.5,  7.5,  12.5,  etc.  The  final  column  gives  the  ap¬ 
proximate  cumulative  distribution  function  defined  on  the  positive  integers 
obtained  from  a  linear  interpolation  between  the  fractiles.  Notice  the  close 
agreement  between  the  exact  and  approximate  cumulative  probabilities  in  the 
tail  of  the  distribution. 

6.  The  First  Four  Moments  of  the  LPG  Distribution 

Knowledge  of  the  moments  of  a  complex  distribution  can  be  utilized  in  a 
variety  of  ways.  The  moments  can  be  used  to  estimate  the  distribution  para¬ 
meters  or  to  approximate  the  distribution  itself.  We  derive  the  first  four 
central  moments  (moments  about  the  mean)  of  the  LPG  distribution. 

The  distribution  of  Z(t),  the  number  of  units  demanded  in  time  t,  is 
negative  binomial  with  parameters  q=0,  p*l-0  and  n*ct.  From  Kendall  and  Stuart 
[7],  the  first  four  cumulants  of  the  negative  binomial  distribution  are  given 
by 

2  3  2  4 

K1  *  nq/p,  K^  *  nq/p  ,  K^  *  nq(l+q)/p  and  -  nq(l+4q+q  )/p  . 

The  first  three  cumulants  are  equal  to  the  first  three  central  moments,  respec¬ 
tively,  while  the  fourth  central  moment,  is  given  by 
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Table  4.  The  Sealed  Negative  Binomial  Approximation  to  the  LPG  Distribution 


(Parameters 

are  the  same 

as  those 

of  Table 

2). 

Exact  Probabilities 

Approximate  Probabilities 

X 

b  (x) 

H(x) 

H(x) 

h  (x) 

Approximate 
Cumul at i ve 

0 

.5000 

.5000 

.1905 

.1905 

1 

.1243 

.624  3 

.1905 

.3810 

2 

.0806 

.  7049 

.6667 

.1905 

.5715 

3 

.0589 

.  7638 

.1174 

.6889 

4 

.0451 

.8089 

.0444 

.  7333 

5 

.0355 

.8444 

.0444 

.  7777 

6 

.0283 

.8727 

.0444 

.8221 

7 

.0228 

.8955 

.8889 

.0444 

.8665 

8 

.0185 

.9140 

5296 

.8961 

9 

.0151 

.9291 

.0148 

.9109 

10 

.0124 

.9415 

CO 

<r 

o 

.9257 

11 

.0101 

.9516 

.01A8 

.9405 

12 

.0083 

.9599 

.9630 

.0148 

.9553 

13 

.0069 

.9668 

.0099 

.9652 

14 

.0057 

.9725 

.0049 

.9701 

15 

.0047 

.9775 

.0049 

.9750 

16 

.0039 

.9811 

.0049 

.9799 

17 

.0032 

.9843 

.9877 

.0049 

.9848 

18 

.0027 

.9869 

.0033 

.9881 

19 

.0022 

.9892 

.0016 

.9897 

20 


0018 


9910 


0016 


9913 


4 


i 

i 

M4  =  K4  +  3K2. 

i 

Hence,  the  first  four  central  moments  (f.f.m.)  of  Z(t),  say  ,  l£i£4,  are 

(26)  y  =  ctO/ (1-0) 

*  1 

(27)  p2  =  ct0/(l-0)2 

(28)  (i  =  ct0(i+0)/(l-0)3 

*  (29)  m4  =  ctO(l+4l>K)2+3ctC)/(l-0)4 

In  order  to  derive  the  f.f.ra.  of  the  LPG  distribution,  we  use  the  following 
relationships  which  can  be  found  in  Parzen  [8],  p.  55:  Let  X  and  Y  be  two 
(dependent)  random  variables.  Then 

(30)  E(Y)  =  E[E(Y| X) ] 

(31)  VAR(Y)  =  E[VAR(YlX) ]  +  VAR[E(Y|X)] 

(32)  p3(Y)  =  E[U3(Y|X)]  +  u3[E(YIX)] 

(33)  p4(Y)  =  E[u4(Y|X)]  +  6E[ VAR(Y|  X) ]  •  VAR[E(Y|X)] 

+  m4[E(Y|X)] 

where 

(34)  m3(Y)  =  E[(Y-E(Y))3] 

(35)  M4(Y)  =  E[[Y-E(Y)J4] 

In  the  context  of  our  problem,  we  interpret  Y  as  Z(t)  and  X  as  t.  It 

< 

follows  that 

J  E(Z(t) )  -  E[E[ Z(t) | t] ] 

,  -  E[ct0/ (1-0) ] 

I 

i 

( 

! 

I 

i 

I 

I 


(36) 


ca0/8(l-0) 


Similarly, 


VAR(Z ( i ) )  =  E[VAR(Z(l)|iJ  =  VAR [E(Z(t)|  t) ] 

=  E[ctO/ ( 1 —0 ) 2 ]  +  VAR[c  r0/ (1-0) 1 
=  coiO/ (ft(  1-0) 2 )  +  (c9)2  a/  (g2  (1-G)2) 

(37)  =  eaO/ ((.->(  1-0)  ]2  {B+eO}. 

Following  the  same  kinds  of  arguments,  one  eventually  obtains 

(38)  H  (Z(O)  =  cap  .  {e2(l+0)  +  2c202} 

*  [3<i-o)r 

(39)  n.(Z(i))  =  caO  .  (B3(l+4G>+02)  +  62cO(3a+l)  +  6Bc2Q2a  + 

4  TbTi-o)14 
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